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Abstract. This paper proposes a new preconditioning scheme for a linear sys- 
tem with a saddle-point structure arising from a hybrid approximation scheme 
on the sphere, an approximation scheme that combines (local) spherical radial 
basis functions and (global) spherical polynomials. Making use of a recently 
derived inf-sup condition 13 and the Brezzi stability and convergence the- 
orem for this approximation scheme, we show that the linear system can be 
optimally preconditioned with a suitable block-diagonal preconditioner. Nu- 
merical experiments with a non-uniform distribution of data points support 
the theoretical conclusions. 



1. Introduction 

Amongst approaches for scattered data approximation on the sphere, the hybrid 
interpolation scheme of von Golitschek & Light [6] and Sloan & Sommariva [12] . 
which employs both radial basis functions and spherical polynomials, seems an 
attractive method, especially when the data is concentrated in some regions (such 
as over mountain ranges and trenches on the Earth's surface), yet relatively sparse 
in other regions. The underlying idea is that radial basis functions can give good 
approximation for rapidly varying data over short distances, whereas the polynomial 
component can more effectively represent smooth variations on a global scale. The 
radial basis functions are centered at data points which are supposed given, and 
the linear combination of radial basis functions is constrained to be orthogonal, in 
a natural sense, to the finite dimensional space of polynomials. 

However, the hybrid scheme poses difficulties in implementation, compared with 
a pure radial basis function approximation, when the number of centers is large. 
In the case of a pure radial basis function approximation with a (strictly) positive 
definite kernel, the resulting linear system has a matrix that is positive definite, 
allowing an iterative solution by the conjugate gradient method, and precondition- 
ing by, for example, the additive Schwarz method, see [7]. For the hybrid scheme, 
in contrast, the linear equations for the relevant expansion coefficients have the 
saddle-point form, see [12], 
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where A € M. NxN is a positive definite matrix arising from the radial basis func- 
tion part of the function approximation, and Q £ Wl NxM is a matrix of spherical 
harmonics evaluated at the data centers, with M < N. (The matrices A and Q are 
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defined properly in Section 2). The saddle-point structure means that the overall 
matrix is not positive definite, and that the conjugate gradient method is no longer 
a suitable iterative solver. More fundamentally, because the matrix has both pos- 
itive and negative eigenvalues, the problem of constructing a good preconditioner 
becomes more delicate. For a thorough review of strategies and challenges for the 
numerical solution of saddle-point problems, see [21 H]. 

In this manuscript we concentrate on the stability of the saddle-point formulation 
of this hybrid scheme, and devise and validate a rapid preconditioned iterative 
solution method for the solution of the equations from the approximation. We 
make use of the Brezzi stability and convergence theorem well known in the context 
of mixed finite elements, along with the new inf-sup condition of |13j to establish 
convergence of the approximation scheme; and then use the inf-sup condition to 
obtain an optimal preconditioner. 

A leading contender amongst the solution methods for equations with the struc- 
ture (JTJ) (see [2]) is the block preconditioning method of [8] employing an approxi- 
mation to the Schur complement 



S := Q T A- 1 Q. 



The use of Schur complement approximations in preconditioners is by now well 
established (see for example [5]) in the setting of mixed finite element methods. In 
particular, Verfurth |14j showed that for the mixed finite element approximation 
of the Stokes flow problem the Schur complement is spectrally equivalent to the 
identity operator (or to the mass matrix or L2 projection matrix in the finite element 
setting), making this a suitable approximation to the Schur complement. Verfiirth's 
proof makes essential use of the Babuska-Brezzi or inf-sup condition (Assumption 
2.1 in [14]), see 0|4]. 

In the present setting, the Schur complement turns out to be spectrally equivalent 
not to the identity operator/matrix, but rather to a specific diagonal but non- 
constant matrix. This spectral equivalence, a main result of the paper, is stated in 
Theorem |U The key ingredient here, in analogy with the known inf-sup condition 
for the Stokes flow problem, is the inf-sup stability condition recently established 
by Sloan and Wendland |13) for the hybrid approximation problem. 

An approximate solver for the primal operator (the radial basis function interpo- 
lation matrix in this case) is also required. This could be provided, for example, by 
the domain decomposition method of Le Gia, Sloan and Tran [7] , or by any other 
preconditioner for the pure radial basis function problem. The resulting block di- 
agonal preconditioner is symmetric and positive definite, hence the preconditioned 
MINRES method ([10], [5]) is applicable to the full problem (which is symmetric 
but not positive definite). 

In Section 2 we formulate the hybrid approximation scheme and establish no- 
tation. Then in Section 3 we describe the inf-sup condition of |13j . and use the 
Brezzi theorem to establish stability and convergence of the scheme. In Section 4 
we turn to preconditioning, establishing there the main spectral equivalence result. 
In Section 5 numerical results are presented (using the primal preconditioner of [7]) 
and we conclude in Section 6. 
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2. Problem formulation 

Let X — Xn = {xi, X2, . . . ,xjv} be a set of N distinct points on the sphere § d 
in K d+1 . Using these as centers we define a radial basis function approximation 
space 

N 

X Xn = Xn ■= {^a 4 0(-,Xi) : at, . . . , a N G R} 

i=l 

with a suitable kernel function <f>. The kernel is assumed to be (strictly) positive 
definite, that is 

N N 

^^a i <£(x i ,x i )a j > 
z=l 3=1 

for every set of points Xn = {xi, . . . ,xjy} € § d and for all N e N, with equality 
for distinct points Xj only if a\ = ot2 — ■ ■ . = O/v =0. 

The native space A/^ is defined as the completion under the inner product 

(2) (^ai^-.x^^a^-jXj-))^ = ^^o^^Xi,^-) 

» 3 » 

of the linear space 

N 

F :={^a^(-,x J ),a,-GM,x,-GS d ,i = l,...,Ar,iVeN}, 
where we insist that the points Xj are distinct. The norm is as usual defined by 

II • IU = (•' ')0 /2 - 

It is well-known that A^ is a reproducing kernel Hilbert space (see [T]) with the 
reproducing kernel <p(-, •). That is 

0(x,y) = 0(y,x), x,yeS d 
0(-,y) G AA , ye§ d 

and for / e A^ 

(3) (/,<K-,y)>0 = /(y), yes d - 

The kernel function needs to be positive definite in order that the inner product 
(•, satisfy the positivity axiom for an inner product. Equivalently, the matrices 
Ax defined by 

(4) (Ax)^ := <Kxi,Xi), i,j = l,...,N 

are positive definite as well as symmetric for every X and every N £ N. 

Taking now a fixed N € N and a fixed set Xn C S d , we may define the usual 
radial basis function interpolant to a continuous function / on S d by 

N 

/jv(x) = ^aj0(x,Xj), 

3=1 

where ai , . . . , a at are such that 

/iv(xi) = /(xj), i = l,...,iV, 
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which we may write as 

N 

^2<f>(xi,Xj)aj = f(xi) for i = 1,...,N. 
3=1 

That is, the vector a — (ai, . . . , «at) t of coefficients satisfies 

(5) A x a = f x , 
where 

(6) fx :=(/(x 1 ),...,/(x JV )) T . 

The hybrid approximation scheme of von Golitschek & Light [6 and Sloan & 
Sommariva, see [13] . employs not only the radial basis functions, but also spherical 
polynomials of total degree up to some conveniently chosen L > 0. We define 

V L = span{Y A * : k = 1, . . . , M(d, £), £ — 0, . . . , L}, 

where Y^k is a spherical harmonic of degree £, that is, the restriction to S d of a 
homogeneous harmonic polynomial in M. d+1 of degree £, and M(d, £) is the dimension 
of the space spanned by the spherical harmonics of degree £. Then Vl is the 
set of spherical polynomials of degree < L. We shall assume that {Y^f. : k = 
1 , . . . , M(d,£), £ = 0, f , . . .} is an orthonormal set with respect to the usual L2 
inner product, that is 

y«,fc(x)y«',fc'(x)da;(x) = S e ^S k ^,, 

where doj(x) denotes surface measure on S d . Then it is well known that {Y^jt : k — 
1, . . . ,M(d,£), £ = 0, 1, . . .} is a complete orthonormal basis for L2(S> d ). 
For a given function /, the hybrid approximation scheme is then to find 

N 



(7) 


ujv,i(x) = y"^aj<ft(x,Xj) E X N 
3=1 


and 






L M{d,l) 


(8) 






1=0 k=l 


such that 




(9) 


WjV,i(Xi) +PN,L( x i) = /( x i)i 



or equivalently, 

N L M{d,i) 

(10) ^a^xj,xj) + ^ ^ ^.fc^feCxi) = /( Xi ), i=l,...,N, 
j=i e=o k=i 

which is to be solved subject to the side condition 

N 

(11) 5^a ig (x,-) = VgeT'i. 

The condition (fTTj) is equivalent, via ([3]), to {q,UN,L)<t> — for all g € Pl, forcing 
the radial basis function component to be A/^-orthogonal to Vl- It also ensures 
that the defining linear system is square and symmetric. 
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The conditions ([T0|) . (fTT j) can also be seen to be those which derive from the 
solution of the constrained optimization problem 

\\\un,l - /||| subject to (Ye, k ,u N , L )ct> = 



mm } 



for all £ = 0, . . . , L, k — 1, . . . , M(d,£), the coefficients /3^fc being the Lagrange 
multipliers in the Lagrangian 

L M{d,l) 



£ = 2 W UN - L _ /III + z)2 z)2 Pe,k(Ye,k,u N , L )<p 

t=0 fe=l 



= — ||UJV,Z, - f\\% + (jPN,L,UN,L)tj,. 

This is therefore another way of expressing the hybrid approximation problem. 

By choosing the spherical harmonic functions as the basis for Vl in (jlljl . we can 
write (flT)f . ([TT j) as a so-called 'saddle-point' linear system of equations 



(12) 



Ax 

XX 



}X,L 






a 




' fx ' 




. P _ 








where Ax is defined by (|4]), a is the vector of coefficients a.j as defined above, j3 
is a vector containing the coefficients fy^ for k = 1, . . . , M(d, £), £ = 0, . . . , L, and 
Qx,l is the N x M matrix defined by 

(13) {Qx, L )i,tk :=Y t , k (xi), i=l,...,N, k = l,...,M(d,£), £ = 0,...,L, 
and 

M :=Y^M{d,l) = dim (V L ). 
1=0 

In the present application we need to prescribe more precisely the nature of the 
kernel 0(x, y). In the first place we shall assume that it is zonal, meaning that 

<Xx,y) = $(x-y) 

for some function $ e C[— 1, 1], where x • y denotes the Euclidean inner product in 
K rf+1 . More precisely, we shall assume that </>(x, y) has an expansion of the form 

00 M(d,e) 

(14) 0(x,y)=X] a e Y etk (x)Y eik (y), 

e=o k=i 

with a 1 > for all £ > 0. That the expansion is zonal follows from the addition 
theorem for spherical harmonics, 

M(d,l) 

X; ^(x)^ fc (y) = ^^PKd + l,x-y), 



fe=i 



where P^ (d + 1, z) is the Legendre polynomial of degree i in dimension d + 1 nor- 
malized to p? (d + 1, 1) = I, and w,i is the total surface measure of S d , 



duj(x). 



In this situation it is well known that the inner product in N$ can be written as 

00 M(d,l) _ _ 

Ut.kVt.k 



(15) 



=0 fe=l 
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where 

ue„k = / u(x)ye,fc(x)dw(x). 
Indeed, as a special case of (|T5|) we find 

oo M(d,l) - oo M(d,<) 

</,*(-,y)>, = £ E ^^ = E £ /^,*(y) = /(y), 

£=0 fc=l °^ 1=0 fc=l 

thus verifying the reproducing kernel property 
If we further assume that for large £ 

(16) a*~(£+l)- 2s , 

then it follows from (|15l) and (|T6|) that the native space A/^ is equivalent to the 
Sobolev space H s (S d ) with inner product 

oo M(d,l) 

(17) =J2 £ (^+l) 2s ^,fc^,fc- 

fcO fc=l 

Technicalities aside, we remark that the essential difficulty in analysing the hy- 
brid approximation is that Xjy and Vl are both subsets of 7V^, and that in an 
appropriate sense both can approximate M$ as N or L tend to oo. We need the inf- 
sup condition, now to be introduced, to allow both subspaces to coexist comfortably 
within the one approximation. 

3. INF-SUP CONDITION, AND THE BREZZI THEOREM 

Typically, uniqueness of the solution and optimal error estimates for saddle- 
point problems follow from so-called inf-sup conditions together with appropriate 
coercivity of the primal operator. For us an essential tool will be the following inf- 
sup theorem proved in [13]. In this theorem hx, for a given point set X — Xx C E> d , 
is the mesh norm, defined by 

hx '■= sup inf cos (x-xj). 

In words, hx is the maximum geodesic distance from a point on E> d to the nearest 
point of X. 

Theorem 1. Let (f> be a kernel satisfying tlM and 116\) for some s > d/2. There 
exist constants 7 > and t > depending only on d and s such that for all L > 1 
and all Xx = {xi,...,xat} C S d satisfying hx < t/L the following inequality 
holds: 

mf sup j, > 7. 

pev L \{o} vex N \{0} \\ v U\\pU 

To use this to prove stability we start with the following well-known theorem 
from Brezzi [4], which is at the heart of most analyses of mixed finite elements. 

Theorem 2. Let H and J be real Hilbert spaces, 0(^1,^2) a continuous bilinear 
form on H x H and b(ip,£) be a continuous bilinear form on J x H. Let {Hx ■ 
N G N} and {Jl : L G N} be sequences of subspaces of H and J respectively. Set 

K = {neH : b(0, rf) = V0 G J}, K n ,l = {77 G H N : b{6, ??) = \/0 G J L }- 
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If 

(18) 370 > such that a(rj, rj) > 7o|Mlif Vry € K U .STa^l, 



an. 

3 7 i > such that sup > 7i||6»||, 7 V0 e J 

(19) and sup > 7i||6»||j V0 e J L , 

nen N \{o} WvWh 

then for every £\ € -ff' and £ 2 € J' and every N, L > £/ie discrete problem of 
finding £,n,l € iijv V'jv.l G </l swc/i £/ia£ 

a(^N,L,v) + b(ipN,L,V) = (£i,Tj) \/r]EH N 
b{6,Z N , L ) = (£ 2 ,0) WeJ L 
has a unique solution, and there exists a constant C — C(7o,7i) > such that 

U-£n,l\\h + M-1>n,l\\j<c( inf ||£-£ v || ff + inf 

where t; € H and ip € J are defined by 

a(tv) + K^v) = (4,0 Vjjeff, 
6(0,0 = (4,0) We J. 

To apply this theorem, we first observe that the hybrid approximation with 
its defining equations (jlOp and (1111) can be written, using the reproducing kernel 
property ([3]), as the problem of finding ujv,l € Ajv and j?jv.l € such that 

(20) (uN,L,v)<t> + (PN,L,v)<t> = (f>v)<t> VneAV, 

(21) (ff,uw,i>*=0 VgeP L . 

To use Theorem [5] we take i7 = A^, J = Vl, Hn = %n and Jl = Vl, with 
the inner product on A/^ being defined by @, and the bilinear forms a(-, •) and 
&(•, •) both equal to the A/^ inner product. The coercivity condition (fT8t is trivially 
satisfied on the whole space H = with 70 = 1 since 



a(u, u) = (u, u),), 



it 



2 



The existence of a constant 71 independent of N and L satisfying p^|) is ensured 
by Theorem [1] provided /ix < t/L. The first part of Theorem [2] then confirms that 
the solution of the system (f20| and (f2Tj) exists and is unique provided /ix < t/L. 
The last part of that theorem defines the comparison quantities: it defines ul G A^ 
and pl ^ Vl such that 

(uL,v)<f + (pl,v)4> = (f,v)<e V^eA/^, 
(q,u L )4, = VqeT'L- 

The second equation says that u_l is orthogonal to the space Vl- In principle the 
orthogonality is with respect to the A^ inner product, but because of the zonal 
property of the kernel it is easy to see that this is the same as the L 2 orthogonal 
projection. Indeed, from (11 5|) we have 

(q,u L )^ = 0VqeV L =>■ (u £ )*,fc = 0for*e[0,L] (q,u L ) L2 = OVq € Vl- 
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The first of the latter set of equations then becomes, on specialising the choice of 
rj to q e Vl , 

(pL,q)4> = (/,<?> yqeV L , 

thus pl is the orthogonal projection of / on the subspace Vl with respect to cither 
the L2 or the Ma inner products. We write this orthogonal projection as PlJ ■ Now 
we can write pl = -Pl/, and Ul = f — Pl = / — Phf ■ Theorem [2] with £ = ul and 
ip = Pl now gives the following convergence result, recovering a result obtained by 
a direct argument in [13]. Note that even though we have taken J = Vl in the 
theorem, the constants 70 and 71 do not depend on L, and hence neither does C . 

Theorem 3. Let <fi be a kernel satisfying fi^| ) and 116\) for some s > d/2. There 
exist constants C > and r > depending only on d and s such that for all L > 1 
and all X = Xn = {xi, . . . , xjv} C § d satisfying hx < t/L the solutions of ©j© 
and © satisfy 



(22) 



11/ - u N , L - Pn ,l\U < CJnf \\(f - P L f) - ZnU- 



Explicit error bounds in the L2 norm for / 6 H s and / 6 iJ 2s can then be 
obtained as in [IB] . 

4. Preconditioning 

Now we turn our attention to the linear algebra aspects of the hybrid approxi- 
mation described in Section 2. 

We have noted already that the hybrid approximation can be written as the 
linear system ([H]), with A x ,Qx,l and f defined by © , CE3]) , and @. 

The solution of saddle-point linear systems such as (p} has received much atten- 
tion in recent years - see [2 for an overview of possible approaches. In particular, 
it was shown in [8] that a suitable preconditioner for the saddle point system 



(23) 

with positive definite A is 

(24) 

where 



A Q 
Q T 



A 
S 



S = Q T A- 1 Q 

is the Schur complement. This is because of the remarkable fact that the product 
of (|2"3"|) by the inverse of (IM|) is a diagonalisable matrix with just three distinct 
eigenvalues, namely 1, (1 ± ^/b)/2. Thus an appropriate Krylov subspace iteration 
such as MINRES (see [10]) will converge in just three iterations. While this is 
generally not a practical preconditioner, an approximate preconditioner of the form 



(25) 



A 
S 



where A is a preconditioner for the problem ([5]) involving only A, and S is a suitable 
Schur complement approximation, will lead to rapid convergence. 
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In the present work we shall assume that an approximate preconditioner A for 
A = Ax is already available; one example would be the domain decomposition 
preconditioner from [7J - this is the one we employ in the numerical results presented 
in the next section. Our interest here is in finding an appropriate approximation 
to the Schur complement Sx = Qx l-^jc Qx,l- We shall see that this is handed to 
us by the inf-sup result in Theorem [TJ 

That inf-sup condition can be stated as 

SU P ii ii > 7iblU <G V L , 
v£X N \{0} \\v\\<t> 

provided hx < t/L. With the help of the Cauchy-Schwarz inequality, this can be 
strengthened to a two-sided inequality, 

(26) \\pU= sup > sup v > jiWpU VpeVl, 

«eA/^\{0} \\ v \\<t> vex N \{o} \\v\\<t, 

provided hx < t/L. To find an equivalent matrix expression, we write p € Vl and 
v € Xn as 

L M(d,l) N 

P = ^2 ^2 /3e,kYe, k , v = an<f>(; x»). 

£=0 fe=l i=l 

With the help of the reproducing property ([3]) , we find 

L Af(<M) at 

£=0 fc=l i = l 
L M(d,£) N 

i=a k=i i=x 

and 

N N 

\\v\U = (v,vfj 2 = C£J2 a ^ x ^ 1/2 = (oFAxa) 1 ' 2 . 

i=l 3=1 

Also, with the aid of (fT5|) we obtain 

Pi.k 



^=0 fe=l 

where is the M x M diagonal matrix given by 

(27) {ki)tk,i'k' = 8w§kk'/a,e- 
Thus in matrix terms (|26p can be written as 

(28) (/? T A L #> sup f®*' L< * 1 >li(P T A L ^ V/3eM M . 

QeE«\{o} (a Ax a) = 

Because Ajf is symmetric and positive definite, the central term can be simplified 

_x 

by the substitution a = A x 2 a, making it 

sup ^^f^ = {P T Q x , l A-JQx,lP)K 
aGR N \{0} (a J a)a 
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with the last step following because the supremum over a is clearly achieved by 
a = (p T Q x l^x 2 ) T ■ Thus in matrix terms (j2f))) can be expressed as 



(29) p I A L p>p 1 Q 1 XiL Ax i Qx,LP>liP J A L /? V/3 e 

Through the above arguments we have established the following theorem. 

Theorem 4. Let (f> be a kernel satisfying and U6\) /or some s > d/2, and let 
Ax and Qx,L be given by and A13\) . For all L > 1 and hx < t/L, where r is as 
in Theorem[ll the Schur complement Sx = Qx l^x Qx,L is spectrally equivalent 
to the diagonal matrix given by \2ty . 

It follows from the theorem that our practical recommendation for the hybrid 
problem is a preconditioner of the form 



A 
K L 



(30) 

where A is an approximation to A, and is defined by ([27 

5. Numerical examples 
We will use the following kernel 



0(x, y) = V(|x - y|) = V(V 2 - 2x • y)> *,ye§ 2 , 

where the radial basis function ip(r) is one of the three choices 

Mr) = (l~r)l, Mr) = (1 - r)^.(4r + 1), ^ 2 (r) = (1 - r)^(35r 2 + 18r + 3) 

with = x for x > and otherwise. Note that -00 G C°(R 3 ),-0i G C 2 (R 3 ) 

and ip2 € C 4 (K 3 ) and each is positive definite (see [TS])- We comment that we have 
not here employed any scaling of the compactly supported RBFs. Using functions 
with smaller support would improve matrix conditioning, but because it would also 
reduce approximation accuracy we have chosen not to use any scaling. It is shown 
theoretically in [9] and verified numerically in Figure[T]that the coefficients ae in the 
expansion of the kernel <f> defined from ipi are of order (1 +£)~ 5 . This is consistent 
with (|I6[) but because the constants in the equivalence are large we have chosen to 
work directly with the Fourier-Legendre coefficients ae, which being f-dimensional 
integrals are easily evaluated numerically. For the sphere § 2 , we have 

a e = 2nj $(t)P e (3,t)dt, where $(t) = ip(>/2 - 2b). 

For the preconditioner, as described in [7], the matrix A is preconditioned using 
a domain decomposition technique. First, given a fixed parameter < v < tt, an 
appropriate set of centers {pi, . . . , pj} C X is chosen over the whole sphere so that 

mincos _1 (pi • pj) > v. 

Second, with another fixed parameter < fi < it/ 3, we decompose the point set X 
into a collection of smaller sets Xj, for j = 1, . . . , J, defined by 

Xj := {x 6 X : cos _1 (x • pj) < /i}. 

The sets Xj with cardinality to,-, for j = 1, . . . , J, may overlap each other and must 
satisfy Uj =1 Xj = X. The restriction operator from to M. mj is denoted by Rj 
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Figure 1. Numerical values of (£+ l) 5 a t 

and the extension operator from R m ' to is Rj . Given a vector r £ K , the 
action of the preconditioner A is given by 

3=1 

where the matrix Aj is the restriction of the full matrix A on the subdomain Xj 
(see [7] for more details). 

By Theorem [4] the diagonal matrix defined by (|2T[) is spectrally equivalent to 
Sx = Q x .l^x 1 Qx,l- The block diagonal preconditioner for (|2"3")l is therefore the 
matrix ([50)1 . 

The results here are for interpolation of the function 

f(x, y, z) = exp(x + y + z) + [0.01 - x 2 - y 2 - (z - l) 2 ] 2 +1 

consisting of a smooth first term and a second term whose support is a cap of 
Euclidean radius 0.1. 

Using each of the kernel functions (j> obtained from the radial basis functions 
ipm, m = 0,1,2, we employ N = 2000, 4000, 8000, 16000 and 32000 points, and 
maximum polynomial degree L = 0,5, 10, 15, 20, 25. In each thousand of 

the points were generated in a cap about the z axis subtending an angle of 0.1 
radians at the origin and the remaining points distributed outside this cap. The 
Saff-Kuijlars equal area algorithm described in [TT] was used to generate these 
points in the following manner. Firstly the Saff-Kuijlars points are generated on 
the whole sphere and those in the cap region are discarded. Then a similar equal 
area construction only for the cap is used to generate 1000 points in this region. 
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m 


N 


2000 


4000 


8000 


16000 


32000 





L = 





207 ( 14) 


228 ( 62) 


273 ( 297) 


294 ( 1229) 


367 ( 6277) 




L = 


5 


813 ( 57) 


1116 (304) 


1531 (1662) 


1996 ( 8297) 


2535 (43513) 




L = 


10 


1053 ( 77) 


1488 (413) 


2240 (2463) 


3016 (12573) 


4001 (67005) 




L = 


15 


1080 ( 85) 


1492 (498) 


2208 (2475) 


3037 (12788) 


4270 (71728) 




L = 


20 


1126 ( 99) 


1466 (464) 


2032 (2334) 


2777 (17227) 


3832 (66025) 




L = 


25 


1262 (123) 


1467 (499) 


1976 (2382) 


2562 (11069) 


3461 (58451) 


1 


L = 





1041 ( 75) 


790 ( 218) 


695 ( 755) 


912 ( 3954) 


995 (17436) 




L = 


5 


2654 (193) 


4293 (1194) 


3564 (3889) 


3172 (13736) 


2321 (40754) 




L = 


10 


3647 (280) 


4370 (1237) 


4047 (4488) 


5335 (42445) 


3381 (59228) 




L = 


15 


3267 (265) 


4308 (1281) 


4021 (4594) 


3792 (16908) 


2808 (48999) 




L = 


20 


2883 (264) 


2923 ( 935) 


3329 (6906) 


3977 (17987) 


>24 hours 




L = 


25 


2659 (265) 


2837 (1006) 


3295 (4048) 


2756 (12735) 


3461 (61280) 


2 


L = 





1079 ( 82) 


1578 (457) 


1903 ( 2164) 


2743 (12233) 


2671 (48177) 




L = 


5 


1978 (152) 


2594 (757) 


14424 (16476) 


>24 hours 


>40 hours 




L = 


10 


2205 (176) 


3130 (929) 


13569 (15603) 


>24 hours 


>40 hours 




L = 


15 


2402 (205) 


2803 (873) 


11970 (14053) 


>24 hours 


>40 hours 




L = 


20 


1796 (169) 


1869 (624) 


8710 (10687) 


13551 (63016) 


>40 hours 




L = 


25 


1615 (168) 


1758 (629) 


6417 ( 8239) 


10299 (48420) 


>40 hours 



Table 1. MINRES iteration count (CPU time) without preconditioning 



The number of MINRES iterations and the CPU time in seconds required for 
convergence to a residual norm tolerance of 10 -9 are tabulated for the unprecondi- 
tioned case in Tabic Q] and for the preconditioner introduced here in Table [2j The 
computer code is written in Fortran 90, compiled with the Intel compiler and run 
on a single core of an SGI Altix XE320 with two Intel Xeon X5472 CPUs. 

The preconditioning is seen to be effective: as anticipated from the theory above, 
the number of iterations remains approximately constant over all choices of N for 
each degree L. Indeed, for each TV aside from the simple case L = (in which the 
radial basis function approximation matrix is only supplemented by one row and 
one column), the iteration counts grow only slowly with increasing L. 

In order to determine the descriptiveness of the bound (|29|) we have also com- 
puted the generalised eigenvalues A; of the pencil Qxl^Qxl — for the 
case N = 4000 points, for m = 0, 1 and for L = 5, 10, 15, 20, 25. Note that the 
generalised eigenvalues are exactly the eigenvalues of A£ 1 (Q3f l&xQx.l)- The 
minimum and maximum computed eigenvalues are given in Table [3J It is notice- 
able how close the largest eigenvalue is to the analytical upper bound of 1 and that, 
although the lowest eigenvalue does decrease for larger L, it remains reasonably 
close to 1. (Note that for fixed X and increasing L the inf-sup condition must 
eventually break down.) 

6. Conclusions 

By employing a recent inf-sup stability result of Sloan and Wendland, we have 
derived an effective preconditioned iterative solver for the hybrid radial basis func- 
tion and spherical polynomial approximation scheme of Sloan and Sommariva. The 



PRECONDITIONING FOR APPROXIMATION ON THE SPHERE 



13 



m 


N 


2000 


4000 


8000 


16000 


32000 





L = 





31 (6) 


39 (18) 


30 ( 51) 


29 (225) 


39 ( 959) 




L = 


5 


59 (11) 


71 (31) 


58 ( 96) 


62 (465) 


75 (1802) 




L = 


10 


70 (13) 


88 (39) 


70 (116) 


71 (532) 


89 (2120) 




L = 


15 


76 (15) 


93 (43) 


76 (128) 


76 (572) 


96 (2279) 




L = 


20 


83 (17) 


98 (48) 


80 (138) 


82 (780) 


99 (2393) 




L = 


25 


95 (20) 


97 (50) 


80 (142) 


84 (649) 


104 (2505) 


1 


L = 





43 ( 8) 


75 (34) 


35 ( 66) 


29 (227) 


47 (1161) 




L = 


5 


76 (15) 


128 (57) 


83 (138) 


74 (557) 


105 (2559) 




L = 


10 


91 (18) 


148 (67) 


98 (163) 


94 (705) 


136 (3296) 




L = 


15 


98 (19) 


168 (78) 


96 (163) 


100 (758) 


148 (3615) 




L = 


20 


107 (22) 


170 (83) 


103 (180) 


103 (788) 


153 (3754) 




L = 


25 


106 (23) 


174 (89) 


103 (185) 


113 (870) 


161 (3972) 


2 


L = 





64 (13) 


149 ( 61) 


46 ( 81) 


30 ( 243) 


61 (1419) 




L = 


5 


95 (19) 


157 ( 70) 


88 (151) 


103 ( 797) 


140 (3380) 




L = 


10 


95 (19) 


171 ( 71) 


102 (176) 


111 ( 861) 


146 (3425) 




L = 


15 


112 (23) 


187 (123) 


113 (199) 


118 ( 923) 


165 (3614) 




L = 


20 


119 (25) 


197 ( 86) 


115 (207) 


133 (1052) 


196 (4301) 




L = 


25 


125 (28) 


201 ( 91) 


119 (220) 


131 (1046) 


203 (4433) 


Table 2. 


MINRES iteration count(CPU time) with precondition- 


ing A and S 










m 


L 


5 


10 


15 


20 


25 







0.9987434 


0.9899326 


0.9623012 


0.9068357 


0.8348191 




^max 


0.9997653 


0.9997658 


0.9997674 


0.9997753 


0.9998099 


1 


-^min 


0.9999955 


0.9999125 


0.9993989 


0.9973949 


0.9908182 




-^max 


0.9999986 


0.9999986 


0.9999986 


0.9999986 


0.9999989 



Table 3. Extreme eigenvalues of Q^^A^Qxx - AA L for N = 4000 



preconditioner requires only a good approximation for the radial basis function in- 
terpolation problem and a simple diagonal scaling matrix for the Schur complement 
based on the Fourier-Legendre coefficients of the kernel function used in the radial 
basis function interpolant. We have established theoretically that the precondi- 
tioner for the Schur complement is optimal. 

Acknowledgements. The support of the Australian Research Council is grate- 
fully acknowledged. 
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